package org.openscales.proj4as.proj {

	import org.openscales.proj4as.ProjPoint;
	import org.openscales.proj4as.ProjConstants;
	import org.openscales.proj4as.Datum;

	/**
	 <p>NEW ZEALAND MAP GRID projection</p>
	 
	 <p>PURPOSE:	Transforms input longitude and latitude to Easting and
	 Northing for the New Zealand Map Grid projection.  The
	 longitude and latitude must be in radians.  The Easting
	 and Northing values will be returned in meters.</p>
	 
	 <p>ALGORITHM REFERENCES</p>
	 
	 <p>1.  Department of Land and Survey Technical Circular 1973/32
	 http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf</p>
	 
	 <p>2.  OSG Technical Report 4.1
	 http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf</p>
	 
	 <p>IMPLEMENTATION NOTES</p>
	 
	 <p>The two references use different symbols for the calculated values. This
	 implementation uses the variable names similar to the symbols in reference [1].</p>
	 
	 <p>The alogrithm uses different units for delta latitude and delta longitude.
	 The delta latitude is assumed to be in units of seconds of arc x 10^-5.
	 The delta longitude is the usual radians. Look out for these conversions.</p>
	 
	 <p>The algorithm is described using complex arithmetic. There were three
	 options:<br />
	 * find and use a Javascript library for complex arithmetic<br />
	 * write my own complex library<br />
	 * expand the complex arithmetic by hand to simple arithmetic</p>
	 
	 <p>This implementation has expanded the complex multiplication operations
	 into parallel simple arithmetic operations for the real and imaginary parts.
	 The imaginary part is way over to the right of the display; this probably
	 violates every coding standard in the world, but, to me, it makes it much
	 more obvious what is going on.</p>
	 
	 <p>The following complex operations are used:<br />
	 - addition<br />
	 - multiplication<br />
	 - division<br />
	 - complex number raised to integer power<br />
	 - summation</p>
	 
	 <p>A summary of complex arithmetic operations:<br />
	 (from http://en.wikipedia.org/wiki/Complex_arithmetic)<br />
	 addition:       (a + bi) + (c + di) = (a + c) + (b + d)i<br />
	 subtraction:    (a + bi) - (c + di) = (a - c) + (b - d)i<br />
	 multiplication: (a + bi) x (c + di) = (ac - bd) + (bc + ad)i<br />
	 division:       (a + bi) / (c + di) = [(ac + bd)/(cc + dd)] + [(bc - ad)/(cc + dd)]i</p>
	 
	 <p>The algorithm needs to calculate summations of simple and complex numbers. This is
	 implemented using a for-loop, pre-loading the summed value to zero.</p>
	 
	 <p>The algorithm needs to calculate theta^2, theta^3, etc while doing a summation.
	 There are three possible implementations:
	 - use Math.pow in the summation loop - except for complex numbers
	 - precalculate the values before running the loop
	 - calculate theta^n = theta^(n-1) * theta during the loop
	 This implementation uses the third option for both real and complex arithmetic.</p>
	 
	 <p>TEST VECTORS</p>
	 
	 <p>NZMG E, N:         2487100.638      6751049.719     metres<br />
	 NZGD49 long, lat:      172.739194       -34.444066  degrees<br />
	 <br />
	 NZMG E, N:         2486533.395      6077263.661     metres<br />
	 NZGD49 long, lat:      172.723106       -40.512409  degrees<br />
	 <br />
	 NZMG E, N:         2216746.425      5388508.765     metres<br />
	 NZGD49 long, lat:      169.172062       -46.651295  degrees</p>
	 
	 <p>Note that these test vectors convert from NZMG metres to lat/long referenced
	 to NZGD49, not the more usual WGS84. The difference is about 70m N/S and about
	 10m E/W.</p>
	 
	 <p>These test vectors are provided in reference [1]. Many more test
	 vectors are available in
	 http://www.linz.govt.nz/docs/topography/topographicdata/placenamesdatabase/nznamesmar08.zip
	 which is a catalog of names on the 260-series maps.</p>
	 
	 <p>EPSG CODES</p>
	 <br />
	 NZMG     EPSG:27200<br />
	 NZGD49   EPSG:4272<br />
	 <br />
	 http://spatialreference.org/ defines these as<br />
	 Proj4js.defs["EPSG:4272"] = "+proj=longlat +ellps=intl +datum=nzgd49 +no_defs ";<br />
	 Proj4js.defs["EPSG:27200"] = "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +datum=nzgd49 +units=m +no_defs ";
	 
	 <p>LICENSE<br />
	 Copyright: Stephen Irons 2008
	 Released under terms of the LGPL as per: http://www.gnu.org/copyleft/lesser.html</p>
	 
	 **/
	public class ProjNzmg extends AbstractProjProjection {
		private var A:Array=[];
		private var B_re:Array=[];
		private var B_im:Array=[];
		private var C_re:Array=[];
		private var C_im:Array=[];
		private var D:Array=[];
		private var iterations:int=1;

		public function ProjNzmg(data:ProjParams) {
			super(data);
		}

		override public function init():void {
			this.A[1]=+0.6399175073;
			this.A[2]=-0.1358797613;
			this.A[3]=+0.063294409;
			this.A[4]=-0.02526853;
			this.A[5]=+0.0117879;
			this.A[6]=-0.0055161;
			this.A[7]=+0.0026906;
			this.A[8]=-0.001333;
			this.A[9]=+0.00067;
			this.A[10]=-0.00034;

			this.B_re[1]=+0.7557853228;
			this.B_im[1]=0.0;
			this.B_re[2]=+0.249204646;
			this.B_im[2]=+0.003371507;
			this.B_re[3]=-0.001541739;
			this.B_im[3]=+0.041058560;
			this.B_re[4]=-0.10162907;
			this.B_im[4]=+0.01727609;
			this.B_re[5]=-0.26623489;
			this.B_im[5]=-0.36249218;
			this.B_re[6]=-0.6870983;
			this.B_im[6]=-1.1651967;

			this.C_re[1]=+1.3231270439;
			this.C_im[1]=0.0;
			this.C_re[2]=-0.577245789;
			this.C_im[2]=-0.007809598;
			this.C_re[3]=+0.508307513;
			this.C_im[3]=-0.112208952;
			this.C_re[4]=-0.15094762;
			this.C_im[4]=+0.18200602;
			this.C_re[5]=+1.01418179;
			this.C_im[5]=+1.64497696;
			this.C_re[6]=+1.9660549;
			this.C_im[6]=+2.5127645;

			this.D[1]=+1.5627014243;
			this.D[2]=+0.5185406398;
			this.D[3]=-0.03333098;
			this.D[4]=-0.1052906;
			this.D[5]=-0.0368594;
			this.D[6]=+0.007317;
			this.D[7]=+0.01220;
			this.D[8]=+0.00394;
			this.D[9]=-0.0013;
		}

		/**
		   New Zealand Map Grid Forward  - long/lat to x/y
		   long/lat in radians
		 */
		override public function forward(p:ProjPoint):ProjPoint {
			var lon:Number=p.x;
			var lat:Number=p.y;

			var delta_lat:Number=lat - this.lat0;
			var delta_lon:Number=lon - this.long0;

			// 1. Calculate d_phi and d_psi    ...                          // and d_lambda
			// For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians.
			var d_phi:Number=delta_lat / ProjConstants.SEC_TO_RAD * 1E-5;
			var d_lambda:Number=delta_lon;
			var d_phi_n:Number=1; // d_phi^0

			var d_psi:Number=0;
			for (var n:int=1; n <= 10; n++) {
				d_phi_n=d_phi_n * d_phi;
				d_psi=d_psi + this.A[n] * d_phi_n;
			}

			// 2. Calculate theta
			var th_re:Number=d_psi;
			var th_im:Number=d_lambda;

			// 3. Calculate z
			var th_n_re:Number=1;
			var th_n_im:Number=0; // theta^0
			var th_n_re1:Number;
			var th_n_im1:Number;

			var z_re:Number=0;
			var z_im:Number=0;
			for (n=1; n <= 6; n++) {
				th_n_re1=th_n_re * th_re - th_n_im * th_im;
				th_n_im1=th_n_im * th_re + th_n_re * th_im;
				th_n_re=th_n_re1;
				th_n_im=th_n_im1;
				z_re=z_re + this.B_re[n] * th_n_re - this.B_im[n] * th_n_im;
				z_im=z_im + this.B_im[n] * th_n_re + this.B_re[n] * th_n_im;
			}

			// 4. Calculate easting and northing
			var x:Number=(z_im * this.a) + this.x0;
			var y:Number=(z_re * this.a) + this.y0;

			p.x=x;
			p.y=y;

			return p;
		}


		/**
		   New Zealand Map Grid Inverse  -  x/y to long/lat
		 */
		override public function inverse(p:ProjPoint):ProjPoint {
			var x:Number=p.x;
			var y:Number=p.y;

			var delta_x:Number=x - this.x0;
			var delta_y:Number=y - this.y0;

			// 1. Calculate z
			var z_re:Number=delta_y / this.a;
			var z_im:Number=delta_x / this.a;

			// 2a. Calculate theta - first approximation gives km accuracy
			var z_n_re:Number=1;
			var z_n_im:Number=0; // z^0
			var z_n_re1:Number;
			var z_n_im1:Number;

			var th_re:Number=0;
			var th_im:Number=0;
			for (var n:int=1; n <= 6; n++) {
				z_n_re1=z_n_re * z_re - z_n_im * z_im;
				z_n_im1=z_n_im * z_re + z_n_re * z_im;
				z_n_re=z_n_re1;
				z_n_im=z_n_im1;
				th_re=th_re + this.C_re[n] * z_n_re - this.C_im[n] * z_n_im;
				th_im=th_im + this.C_im[n] * z_n_re + this.C_re[n] * z_n_im;
			}

			// 2b. Iterate to refine the accuracy of the calculation
			//        0 iterations gives km accuracy
			//        1 iteration gives m accuracy -- good enough for most mapping applications
			//        2 iterations bives mm accuracy
			for (var i:int=0; i < this.iterations; i++) {
				var th_n_re:Number=th_re;
				var th_n_im:Number=th_im;
				var th_n_re1:Number;
				var th_n_im1:Number;

				var num_re:Number=z_re;
				var num_im:Number=z_im;
				for (n=2; n <= 6; n++) {
					th_n_re1=th_n_re * th_re - th_n_im * th_im;
					th_n_im1=th_n_im * th_re + th_n_re * th_im;
					th_n_re=th_n_re1;
					th_n_im=th_n_im1;
					num_re=num_re + (n - 1) * (this.B_re[n] * th_n_re - this.B_im[n] * th_n_im);
					num_im=num_im + (n - 1) * (this.B_im[n] * th_n_re + this.B_re[n] * th_n_im);
				}

				th_n_re=1;
				th_n_im=0;
				var den_re:Number=this.B_re[1];
				var den_im:Number=this.B_im[1];
				for (n=2; n <= 6; n++) {
					th_n_re1=th_n_re * th_re - th_n_im * th_im;
					th_n_im1=th_n_im * th_re + th_n_re * th_im;
					th_n_re=th_n_re1;
					th_n_im=th_n_im1;
					den_re=den_re + n * (this.B_re[n] * th_n_re - this.B_im[n] * th_n_im);
					den_im=den_im + n * (this.B_im[n] * th_n_re + this.B_re[n] * th_n_im);
				}

				// Complex division
				var den2:Number=den_re * den_re + den_im * den_im;
				th_re=(num_re * den_re + num_im * den_im) / den2;
				th_im=(num_im * den_re - num_re * den_im) / den2;
			}

			// 3. Calculate d_phi              ...                                    // and d_lambda
			var d_psi:Number=th_re;
			var d_lambda:Number=th_im;
			var d_psi_n:Number=1; // d_psi^0

			var d_phi:Number=0;
			for (n=1; n <= 9; n++) {
				d_psi_n=d_psi_n * d_psi;
				d_phi=d_phi + this.D[n] * d_psi_n;
			}

			// 4. Calculate latitude and longitude
			// d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians.
			var lat:Number=this.lat0 + (d_phi * ProjConstants.SEC_TO_RAD * 1E5);
			var lon:Number=this.long0 + d_lambda;

			p.x=lon;
			p.y=lat;

			return p;
		}

	}
}